perf(engine): sparse dynamic solvers — time history and harmonic - #180

…gration

The buckling paths built K_g as a dense n×n matrix and then converted it
to CSC with from_dense_symmetric for the sparse Lanczos eigensolver — the
n² allocation the sparse path exists to avoid. The constraint paths
additionally densified K_ff and ran dense Jacobi.

- geometric_stiffness: build_kg_from_forces_2d/3d and the five shell
  add_*_geometric_stiffness_3d now scatter into a KgTriplets sink that
  keeps only the free×free block and converts to lower-triangle CSC
  (duplicates summed). Element loops emit each unordered local pair once
  (lower triangle) — emitting the full symmetric square would
  double-count off-diagonals through from_triplets' duplicate summation
  (caught by the Euler buckling goldens during development).
- buckling 2D/3D: no dense K_g, no extract/negate roundtrip; the
  constraint path reduces K and -Kg with the sparse triple product and
  uses the same sparse shift-invert Lanczos as the unconstrained path
  (dense Jacobi gone from buckling).

Full suite green (7101), including the buckling goldens and the
sparse-vs-Jacobi parity test (whose manual reference path now densifies
the CSC triplets instead of the other way around).

Stacked on perf/sparse-constraint-transform (#177).

numeric_cholesky used to call permute_symmetric on every factorization:
tripletize + sort_unstable O(nnz log nnz) plus six nnz-sized temporaries,
all of it purely symbolic work. The symbolic phase now stores the
permuted structure (pa_col_ptr/pa_row_idx) and pa_src[p] = index into the
original values array, so a numeric factorization permutes values with a
single O(nnz) gather and no allocations beyond the output.

The pattern-reuse contract is unchanged (callers already had to pass the
same-structure matrix); the doc comment now states it explicitly.

Measured on the 256x256 grid (n=65536, nnz_L=2.02M), 5 repeated numeric
factorizations with the symbolic reused (the P-Delta pattern): ~116 ms ->
~100 ms per factorization (medians of 3 runs; noisy box). The win is
bigger in allocation pressure than in wall time: no per-call triplet
buffers.

Stacked on perf/amd-quotient-graph (#176).

2D paths dispatch at SPARSE_THRESHOLD: small models keep the original
dense code, large models assemble K (and M for Craig-Bampton) as CSC,
reduce constraints as sparse triple products, keep K_II as a sparse
principal submatrix, and factor it once with sparse Cholesky (dense LU
fallback capped by MAX_DENSE_FALLBACK_DOFS). Craig-Bampton interior
modes come from sparse shift-invert Lanczos. Guyan reactions use the
sparse full-K matvec instead of dense K_rf/K_rr extraction.

3D paths stop densifying when constraints are present
(reduce_matrix_sparse instead of reduce_matrix(to_dense)) and get the
same sparse K_II factorization; Craig-Bampton 3D drops the dense O(n^2)
mass assembly.

Blocks involving only boundary DOFs stay dense (nb is small by design);
the small-model 2D path is byte-identical to the old code.

Entries that are analytically zero pick up rounding-level asymmetry from
the different operation order of the sparse products; comparing the
difference against the entry pair itself fails on near-zero couplings.
Scale against max|K| instead.

Above SPARSE_THRESHOLD the 2D/3D co-rotational Newton loops now convert
the extracted free tangent block to CSC, reduce constraints with
reduce_matrix_sparse, and factor with sparse Cholesky — the symbolic
factorization is computed once per solve call and reused across
iterations (pattern-fingerprinted; rebuilt only if the CSC structure
changes). Modified NR caches the sparse factor per increment. Non-SPD
tangents fall back to dense LU, capped by MAX_DENSE_FALLBACK_DOFS.

Below the threshold the dense path is unchanged. New validation test:
30-element cantilever (nf = 90) converges through the sparse path, tip
deflection within 1% of PL^3/3EI, dense/sparse and full/modified NR
parity within 1e-4.

…ment control

FactoredTangent gains a Sparse variant (numeric Cholesky over a cached
symbolic factorization, fingerprinted by the CSC pattern and reused
across all steps/iterations of a solve call). Above SPARSE_THRESHOLD the
constraint reduction moves to reduce_matrix_sparse; non-SPD tangents
fall back to dense LU capped by MAX_DENSE_FALLBACK_DOFS. Factor-once /
solve-many is preserved (arc-length corrector still solves 2 RHS per
factorization). Below the threshold the dense path is byte-identical.

New validation tests: 30-element cantilever (nf=90, sparse) vs
21-element mesh of the same cantilever (nf=63, dense) — arc-length and
displacement-control load factor and tip displacement agree within 2%.

…solvers

SparseSymbolicCache / cached_symbolic / tangent_free_sparse /
solve_tangent_sparse move from corotational.rs (with a near-duplicate in
arc_length.rs) into solver/sparse_tangent.rs. The shared
tangent_free_sparse takes the already-extracted free block (the
arc_length signature); corotational's call sites inline the extraction
the old private copy did internally. Pure move, no behavior change.

Both contact entry points (2D/3D) dispatch at SPARSE_THRESHOLD to the
shared sparse_tangent helpers: CSC conversion, sparse constraint
reduction, sparse Cholesky over a fingerprinted symbolic cache. The
active set toggles penalty entries per iteration, which changes the CSC
pattern — the fingerprinted cache rebuilds the symbolic exactly then and
reuses it while the active set is stable. Below the threshold the dense
path is byte-identical, error messages included.

New parity test: cantilever pressed against a rigid wall through a gap
element, coarse mesh (dense) vs fine mesh (sparse) — gap displacement
and transmitted force agree within 1e-6.

2D/3D fiber Newton loops dispatch at SPARSE_THRESHOLD to the shared
sparse_tangent helpers; modified NR caches the sparse numeric factor of
the increment's first tangent (the analogue of the dense cached_l).
Dense path below the threshold is unchanged. New parity test: meshed
elastic cantilever below/above the threshold, tip deflection vs PL^3/3EI
and dense/sparse/modified-NR parity within 1e-4.

The Ernst-iteration solve (2D and 3D) dispatches at SPARSE_THRESHOLD to
the shared sparse_tangent helpers with the symbolic Cholesky cached
across iterations. New parity test: harp-stayed deck straddling the
threshold, cable tension and deck midspan deflection agree within 2%.

The secant-stiffness iteration (2D and 3D) dispatches at
SPARSE_THRESHOLD to the shared sparse_tangent helpers. New parity test:
soft-clay lateral pile at two mesh densities straddling the threshold
(head deflection agrees within 7%; the secant iteration's mesh
sensitivity is pre-existing and solver-independent — verified by running
the finer mesh through the dense path).

The per-stage solves (2D, 3D, and the staged cable-tension loop)
dispatch at SPARSE_THRESHOLD to the shared sparse_tangent helpers. The
active-element set changes the sparsity pattern between stages; the
fingerprinted symbolic cache rebuilds exactly then. New test: two-stage
construction matches a linear reference solve on meshes both sides of
the threshold, within 1e-6.

Above SPARSE_THRESHOLD the 2D/3D Newton loops assemble the tangent
directly as lower-triangle COO triplets instead of a dense n*n matrix:
the element assemblers (truss/frame corotational, springs) are now
generic over a scatter closure, so the dense path keeps byte-identical
arithmetic. Inclined-support rotation uses the COO triplet helpers from
assembly.rs. The final k_dummy dense matrix is only assembled when
constraints exist (constraint forces need it); otherwise f_int is
accumulated via the triplet drivers.

assemble_corotational_public keeps its dense signature (arc-length still
consumes it). New test: propped cantilever with an inclined roller
straddling the threshold — kinematic restraint check and dense/sparse
parity within 1e-4.

The fiber/elastic/spring assemblers are generic over a scatter closure
(dense path byte-identical); above SPARSE_THRESHOLD the Newton loops
assemble lower-triangle COO triplets and build the CSC free block
directly — no dense n*n tangent, no extract_submatrix.
tangent_free_sparse_triplets moves to the shared sparse_tangent module
(corotational imports it from there now). The constraint-forces tangent
rebuild stays dense via the dense closure, mirroring corotational's
k_dummy.

… control

The sparse path now assembles the co-rotational tangent directly as
lower-triangle triplets via assemble_corotational_triplets_2d (now
pub(crate)) instead of a dense n*n matrix. Springs go through a scatter
closure (dense path byte-identical). The final constraint-forces block
uses compute_constraint_forces_sparse over the unreduced CSC free block
— no dense K anywhere on the sparse path. Inclined-support behavior is
unchanged from the status quo (arc-length never rotated the tangent).

Above SPARSE_THRESHOLD the base stiffness is assembled once as CSC
(sparse_assembly::assemble_stiffness_sparse_2d — the production path
with inclined transforms on triplets — and assemble_sparse_3d with
full-K for reactions) and per iteration the corrected tangent is rebuilt
from cached base triplets plus the Ernst correction triplets, never
touching a dense n*n matrix. The Ernst factor computation is shared by
both paths (ernst_diff_2d/3d helpers) so corrections are numerically
identical; the documented unrotated-dk inclined-support gap is
preserved as-is. Reactions use k_full.extract_block_dense, constraint
forces use compute_constraint_forces_sparse. Dense path below the
threshold is unchanged.

Above SPARSE_THRESHOLD the base stiffness is assembled once as CSC and
per iteration the tangent is rebuilt from cached base triplets plus the
soil-spring diagonal triplets — no dense n*n clone per iteration.
Constraint forces use compute_constraint_forces_sparse. Dense path
below the threshold is unchanged.

assemble_tangent_stiffness (2D/3D) is now scatter-generic — the dense
path writes through a closure with byte-identical arithmetic, the sparse
path (ns >= SPARSE_THRESHOLD) pushes lower-triangle triplets, rotates
them for inclined supports with the COO helpers, and builds the CSC
free block directly. The solve upgrades to the shared
solve_tangent_sparse with a per-call cached symbolic Cholesky (the old
sparse_cholesky_solve_full recomputed the symbolic every iteration).
The final reaction pass keeps a dense tangent (k_dummy precedent).
New parity tests: cantilever straddling the threshold, plain and with a
45-degree inclined roller.

Above SPARSE_THRESHOLD the contact solver assembles the base stiffness
once as CSC (2D: the production sparse path with inclined transforms on
triplets; 3D: assemble_sparse_3d) and rebuilds each iteration's tangent
from cached base triplets plus the currently-active penalty corrections
(gap blocks, friction outer products, deactivation as negative element
triplets) — a contact pair going inactive simply contributes no triplets
that iteration, so removal needs no explicit subtraction. The
fingerprinted symbolic cache rebuilds only when the active set changes
the pattern. Dense path below the threshold is verbatim.

The branch forked at the #174 merge and cherry-picked the sparse Cholesky
series out of the queue, so it carried those files in their pre-review form
while main gained four rounds of follow-up fixes on top of them. Merging
produced 18 conflict hunks across five files, every one of them a choice
between this branch's older copy and the corrected version on main.

Resolved so that main's corrections survive and this branch's new work does
too:

geometric_stiffness.rs — main. Recovers the warping remap (#179): the
  buggy side indexes a 12x12 kg_global with elem_dofs.len(), which is 14
  when any section declares `cw`, so kg_global[10 * 14 + 4] walks one past
  the end of a 144-element array. Native that is an index panic; on the
  shipped wasm32 build it is an abort across the FFI boundary. Also
  recovers the 1e-30 pruning of the assembled K_g.

buckling.rs — main, wholesale. Checked first that this branch adds nothing
  of its own here: its diff against the fork point is the same sparse-K_g
  cherry-pick main already has, in its earlier form. Main additionally has
  the K-metric Lanczos fix and the curved-shell gate.

sparse_chol.rs — main, wholesale. Recovers pa_fingerprint and the live
  (not debug_assert) reuse guard. Worth stating plainly: this branch's own
  sparse_tangent.rs is the heaviest consumer of that contract, reusing a
  symbolic factorization across every Newton iteration. Taking the other
  side would have removed the guard protecting its own new feature.

arc_length.rs — a genuine merge of both. This branch's sparse dispatch
  (factor_assembled over a cached symbolic) is kept, and so is #174's reuse
  of `residual_s`: the branch shadowed it with a second
  cs.reduce_vector(&residual), which is exactly the recompute #174 removed.
  The final constraint-force block takes this branch's version, which
  dispatches sparse/dense instead of always building the dense final
  tangent.

tests/.../arc_length.rs — both sides. The two added independent fixtures
  and independent tests at the same place; the trailing brace was shared
  context, so this branch's last test needed its own.

Also dropped a `free_idx` binding left unused by the resolution — #174
hoisted it for a dense extract_submatrix that factor_assembled now owns, so
it had no reader left and the lint job would have failed on it.

Full engine suite green in both profiles, 28 test binaries each.

Not fixed here, because it is not this branch's defect
------------------------------------------------------
The review of this PR flagged that the sparse non-SPD fallback in
factor_tangent is unreached by any test. Probing it settled something
larger: NOTHING in the engine suite ever produces an indefinite tangent.

Measured twice, with an unconditional panic! at the head of each branch and
a full suite run:

  sparse non-SPD branch   28 binaries, zero hits
  shared dense LU fallback 4750 tests, zero hits

test_arc_length_detects_limit_point, the one test named for this, never
traces past a limit point. Both its toggle frame and a 12-segment mesh of
the same geometry rise monotonically to lambda = 1.975 and stop:
turned_over=false for both. The reason is in the fixture — it is built from
`frame` elements, which carry moment, so the structure bends rather than
snapping. A von Mises truss of these dimensions would turn over near
P_cr = EA(h/L)^3, about 0.25; the traced path reaches 400 times that
without a limit point. The test cannot report this because its only
snap-through assertion is `assert!(true, ...)` inside an `if`.

So the untested fallback predates this branch: main has the same dense LU
path, equally unreached. This branch adds a sparse route into it and a size
ceiling, neither of which any test touches either.

Fixing it means building a fixture that genuinely snaps — truss elements or
an apex hinge — and giving the existing test an assertion that can fail.
That is a change to main's test suite, not to this PR, and it is better
done where it can be reviewed as what it is.

One correction to the review while I am here: it said the ceiling makes the
post-limit-point path slower than main. That is wrong. Main attempts a
dense Cholesky that fails at O(n^3/3) before its dense LU; this branch fails
a sparse Cholesky instead and pays two O(n^2) conversions. It is faster.
The ceiling only changes behaviour above 6000 DOFs, where main would
attempt an 800 MB dense LU that wasm32 cannot satisfy anyway.

Claude-Session: https://claude.ai/code/session_01L67REmkuj14xQGpVujp5z9

…to HEAD

Above SPARSE_THRESHOLD each stage assembles straight to lower-triangle
triplets over all DOFs (active elements, prestress FEFs, springs,
artificial stiffness) — no dense n*n. One full-K CSC per stage serves
both the free-block solve (tangent_free_sparse_triplets + cached
symbolic) and the K_fr*u_r correction via sparse_cross_block_matvec;
reactions come from a full-K matvec instead of dense K_rf/K_rr
extraction. 3D composes the shared sparse 3D assembler with the staged
artificial-stiffness logic as extra diagonal triplets. The staged cable
loop rebuilds base + Ernst correction triplets per iteration. Dense
path below the threshold is unchanged. New test: staged stayed deck
with cables straddling the threshold (dense/sparse parity within 5%).

…to perf/sparse-time-integration

The feature-gated phase benchmarks referenced renamed struct fields
(SolverNode.z, SolverSupport.dz/dry, SolverNodalLoad.fz/my,
Displacement.uz) and compared an f64 against an integer literal. CI
never compiles this file (manual-bench-phases is off by default), so the
rot went unnoticed.